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ABSTRACT 

Context. With the advent of visible and infrared long-baseline interferometers with more than two telescopes, both the size and the 
completeness of interferometric data sets have significantly increased, allowing images based on models with no a priori assumptions 
to be reconstructed with an aperture synthesis technique. 

Aims. Our main objective is to analyze the multiple parameters of the image reconstruction process with particular attention to the 
regularization term and the study of their behavior in different situations (types of astrophysical objects, telescope array configurations, 
level of noise, etc.). The secondary goal is to derive practical rules for the users. 

Methods. Using the Multi- aperture image Reconstruction Algorithm (MiRA), we performed multiple systematic tests, analyzing 1 1 
regularization terms commonly used. The tests are made on different astrophysical objects, different {u, v) plane coverages and several 
signal-to-noise ratios to determine the minimal configuration needed to reconstruct an image. We establish a methodology and we 
introduce the mean-square errors (MSE) to discuss the results. 

Results. From the ~24000 simulations performed for the benchmarking of image reconstruction with MiRA, we are able to classify 
the different regularizations in the context of the observations. We find typical values of the regularization weight. A minimal (u, v) 
coverage is required to reconstruct an acceptable image, whereas no limits are found for the studied values of the signal-to-noise ratio. 
We also show that super-resolution can be achieved with increasing performance with the (u, v) coverage filling. 
Conclusions. Using image reconstruction with a sufficient (w, v) coverage is shown to be reliable. The choice of the main parameters of 
the reconstruction is tightly constrained. We recommend that efforts to develop interferometric infrastructures should first concentrate 
on the number of telescopes to combine, and secondly on improving the accuracy and sensitivity of the arrays. 

Key words. Instrumentation: interferometers - Techniques: interferometry - Techniques: image processing 



1. Introduction 

Many astrophysical studies require milli- arc second (hereafter 
mas) resolution images at optical wavelengths (visible and in- 
frared), for example, the understanding of the interplay between 
accretion and ejection in the inner part of the disks of young 
stellar objects, the expansion mechanisms in novae just a few 
hours or days after the explosion, and the nature of dust in active 
galactic nuclei. Information at such a high resolution at the op- 
tical wavelengths requires diffraction-limited images with pupil 
sizes of the order of tens to hundreds of meters that can only 
be achieved by interferometrically combining light from sepa- 
rate apertures. The Very Large Telescope Interferometer (VLTI) 
and the Center for High Angular Resolution Array (CHARA) 
are facilities that provide interferometric measurements that can 
be used to reconstruct images of stellar surfaces (e.g., Monnier 
et al. 2007; Haubois et al. 2009; Zhao et al. 2009), binaries (e.g., 
Zhao et al. 2008; Kraus et al. 2009), circumstellar shells around 
evolved stars (e.g., Le Bouquin et al. 2009), and the close envi- 
ronment of young stars (Renard et al. 2010; Kraus et al. 2010). 

Owing to the sparse (u, v) coverage, the image reconstruction 
process is ill posed as there are more unknowns, e.g. the pixels 
of the image, than measurements. The data alone are insufficient 
to reconstruct an unambiguous image and some additional con- 
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straints, so-called the regularizations, are needed to converge to 
a unique and stable solution. Compared to radio interferometry, 
the data are much sparser in optical interferometry; hence, we 
expect the image reconstruction problem to be much more sen- 
sitive to the choice and the tuning of the regularization. Since 
the general study of regularization by Titterington (1985), many 
different methods have been proposed in the literature to adjust 
the regularization level. Since image reconstruction for optical 
interferometry is still in its infancy, it is fundamental to ana- 
lyze the different types of regularization to find those that are 
the most suitable for the different astrophysical problems and to 
be able to tune the weight of the regularization. In this context, 
we carried out systematic tests using the image reconstruction 
algorithm devoted to optical interferometry data developed by 
Thiebaut (2008), called the Multi-aperture image Reconstruction 
Algorithm (MiRA). The analysis of these tests allow us to ex- 
tract some general conclusions and establish practical rules for 
the users. 

The mathematical principles of the image reconstruction 
technique is presented in Sect. 2. The parameters of the simu- 
lated data are presented together with the characteristics of the 
images and the strategy in Sect. 3. The results of the simulations 
are presented and discussed in Sect. 4, with an analysis of the 
role of the different terms and parameters in the image recon- 
struction. Finally, our conclusions are summarized in Sect. 5. 
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2. Principles of image reconstruction from optical 
interferometric data 

We do not intend to provide a formal and precise description of 
image reconstruction in optical interferometry, but instead suffi- 
cient details to ensure that the paper is self-contained. Readers 
interested by a more detailed description of the method should 
refer to Thiebaut & Giovannelli (2010). 

2.1. Data from optical Interferometrlc observations 

The principle of interferometry is to recombine coherently the 
beams from two or more independent telescopes and measure 
the so-called complex visibilities of the fringe patterns produced 
by the interferences. According to the van Cittert-Zernicke the- 
orem, for an ideal interferometer, the complex visibility Vj lt j 2 (t) 
of the fringes produced by the interferences of the telescopes 
j\ and j2 at time t is proportional to the Fourier transform of 
the object brightness distribution I(Vj u j 2 (t)) at spatial frequency 
v jiJi( t ) ~ Bfrjffl/A, where A is the wavelength, and the so- 
called baseline B± - 2 if) represents the separation between the 
two telescopes projected on a plane perpendicular to the line of 
sight (Lawson 2000; Malbet & Perrin 2007). Since the number 
of measurements is finite, to simplify the equations we introduce 
some notation for the m-th measured complex visibility and the 
corresponding spatial frequency, given by 

Vm = V jhm , hm (t m ), (1) 

v m = Bl imj2m (t m )/A, (2) 

where and j2, m are the interfering telescopes and t m is the 
time of observation. 

2.2. Description of the Image model 

The final product of the image reconstruction is an image that 
can be treated as a grid of square pixels. In this context, the ob- 
ject brightness distribution as a function of the position can be 
approximated using the parametrization 

N 

I(0) = Y J x n b n (0), (3) 

n=\ 

where x = {x n }„ =l are the image parameters, e.g. the pixel val- 
ues of the image, and {b n (0)}^ =l is the chosen basis of functions, 
e.g. the response function of each pixel. The image reconstruc- 
tion then consists of estimating the N parameters x that most 
closely fit the interferometric data. In this paper, we chose x n to 
be proportional to the value of the n-th pixel of a sampled image 
and b n (Q) = b(6 - n ), where b(6) is the pixel shape and n the 
position of the n-i\\ pixel; thus 

x n = al(0 n ), (4) 

where a > is a scaling factor such that x is normalized (this 
is required by the interferometric data format, cf. Pauls et al. 
(2005)). With this model, the exact Fourier transform of the 
brightness distribution is given by 

/(v) = J] x n b n (y) = b(v) ^ *n e- l2n6n - v , (5) 

n n 

where b n (v) and b(y) are the Fourier transforms of the basis func- 
tions. In our case, they correspond to the Fourier transform of the 



pixel response function, i.e. the pixel shape. Hence, the model of 
the m-th complex visibility is given by 

L = hVm) = A m , n X n = (A • X) m , (6) 

n 

where A is a matrix with the complex coefficients 

A m , n = b n (v m ) = b(vje- i2 * »- v ™. (7) 

This matrix multiplication performs a linear transformation that 
contains the Fourier transform, the pixel shape, and the sparse 
sampling of the (u, v) plane. 

The main problem in optical interferometry is the small num- 
ber of telescopes (currently up to four or six), which leads to a 
sparse sampling of the spatial frequencies, the so-called (u, v) 
plane (see Fig. 2). Owing to random effects caused by the atmo- 
spheric turbulence, the visibility phase cannot be calibrated and 
the power spectrum and the closure phase are used. This results 
in a partial loss of the Fourier phase information (Thiebaut & 
Giovannelli 2010). 

2.3. Inverting the problem of Interferometrlc Imaging 

Because of the sparse (u, v) coverage and the possible lack of 
other information such as the phase, the reconstruction of an im- 
age obtained by the interferometric data alone is an ill-posed 
inverse problem. It needs additional a priori constraints to be 
recasted into a problem that has a unique and stable solution. 
A general prescription is to express the solution as one that 
minimizes a penalty function / under some strict constraints 
(Thiebaut 2005; Thiebaut & Giovannelli 2010) 

x + = argmin f(x) s.t. x > and ) x n = 1 , (8) 

x ^ n 

with 

fix) = /data W + A* /priorC*) , (9) 

where the so-called likelihood term /dataC*0 measures the dis- 
crepancy between the model and the available data, while the 
so-called regularization term / pr j 0r (x) measures the discrepancy 
with the prior information. In other words, minimizing the like- 
lihood term /dataC*0 enforces the fit with the actual data, while 
minimizing the regularization term / pr j 0r (x) enforces the agree- 
ment with the priors. The so-called hyperparameter [i > is used 
to adjust the relative weight of the constraints set by the mea- 
surements and the ones set by the priors. In Eq. (8), the positiv- 
ity (x > 0) and the normalization Q} n x n = 1) of the brightness 
distribution are also included by default. 

For the image reconstruction, we use MiRA (Thiebaut 2008) 
to find solutions of Eq. (8). MiRA can deal with the various kinds 
of data provided by an optical interferometer and implements a 
number of different regularizations (Thiebaut 2008; Thiebaut & 
Giovannelli 2010). 

2.3.1 . The likelihood term / data and the data model 

We focus on the choice of the priors and the tuning of their pa- 
rameters. It is not within the scope of the paper to deal with 
global optimization issues and the search for a global minimum 
of the penalty function in Eq. (9). We therefore assume that the 
available measurements consist of complex visibilities, i.e. am- 
plitude and phase, in order to have a convex likelihood term 
/dataC*0- If the regularization term / pr i r(^) is also convex, the 
global penalty function f(x) will be convex, which ensures that 
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the solution of Eq. (8) is unique. Current optical interferome- 
ters only provide phase closures and power- spectrum data (i.e. 
the phase of the bispectrum and the squared amplitude of the 
complex visibilities), this means that our assumption will give 
somewhat optimistic results because some Fourier phase infor- 
mation is missing and because the likelihood term /dataW is 
non-convex when dealing with real data. However, as the num- 
ber of simultaneous interfering telescopes increases, the number 
of missing phases becomes much less important and they can 
be reliably derived using self-calibration (Pearson & Readhead 
1984) to achieve a situation similar to the case studied in our 
simulations. Moreover, new interferometers will make use of a 
phase reference source to directly measure the phase of the com- 
plex visibilities (Delplancke et al. 2003). 

However, the OI-FITS standard (Pauls et al. 2005) imposes 
the use of complex visibilities in their polar representation with 
independent error bars. We therefore simulate each measured 
complex visibility as 

Qm = \V m \ + Sg m , and (f m = arg V m + 6(p m , (10) 

where g m and (f m are the measured amplitude and phase of the m- 
th measure, V m the corresponding complex visibility computed 
from the true object brightness distribution, and 6g m and S(f m are 
additive noise terms. In our simulations, the noise terms have 
independent Gaussian statistics such that 

Var(6g m ) = (g m ) 2 Var(^ m ), (11) 

where (g m ) = \V m \is the expected value of the amplitude that can 
be computed from the complex visibility of the true image. This 
particular choice follows approximately the model of Goodman 
(1985). 

In this context, to define the likelihood term we use the local 
approximation (Meimon et al. 2005) 




where r// fjn (x) and r^ m (x) are the two components of the com- 
plex residuals, r m (x) = g m t lcp,n - (A • x) m , respectively, along 
and orthogonally to the measured complex visibility. Given the 
error bars cr g ^ m and o~^ m of the amplitude and the phase of the 
complex visibility, the standard deviations in the components of 
the complex residuals are (Pauls et al. 2005) 

o-//, m = o~ g , m and cr^ m = g m o~^ m . (13) 
2.3.2. The regularization term f 9r \ or 

In our simulations, we test 11 different regularization terms that 
are commonly used in image reconstruction methods and are im- 
plemented in MiRA (see Appendix A for detailed expressions). 

1. Quadratic smoothness, which most closely describes a 
smooth image and helps us to avoid unmeasured high fre- 
quencies. 

2-3. Compactness, which describes compactness in the im- 
age plane and hence smoothness in the Fourier plane (Le 
Besnerais et al. 2008). Two different cases were studied in 
the simulations, with penalties of the second and third orders 
with respect to the distance of the center of the field of view. 



4. Total variation (hereafter TV), which minimizes the total 
gradient of the image and helps us to describe uniform areas 
in the sought image with steep but localized changes (Strong 
& Chan 2003). 

5. £\ smoothness, which is useful for an extended object with 
sharp edges since it is linear for strong gradients. 

6-8. ^-norm with p - 1.5, p = 2 and p = 3. For p > 1, the 
^p-norm regularization tends to produce a smooth image as 
it reduces the variance in the pixels. 

9-11. Maximum entropy methods (hereafter MEM) aims to 
identify the least informative image consistent with the data 
(Gull & Skilling 1984; Narayan & Nityananda 1986). We 
try three types of entropies MEM-sqrt, MEM-log, and MEM- 
prior, respectively. The first two tend to reproduce an image 
with a flux spread across a minimum number of pixels. The 
last one is minimum when the image is as close as possible 
to a prior image. This prior image is the Gaussian that most 
closely reproduces the amplitude visibility data. 

The different regularization terms are expected to behave as 
follows. The positivity and the normalization imposed in all the 
reconstructions are an l\ norm and lead to the sparsity of the 
solution, i.e. to a minimum number of bright pixels to explain 
the data. As most astrophysical images are smooth and compact, 
we expect that the regularizations that can describe these images 
will behave well, i.e. smoothness (quadratic or compactness, 
and TV. The ^-norm regularization with p = 2 (and by exten- 
sion for p > 1 as the regularization has the same behavior) has 
the tendency to force to zero the spatial frequencies that have 
not been measured (according to the Parseval theorem). Since 
the regularization has to interpolate correctly between the data, 
which is closer to the reality, it is not expected to give good re- 
sults. Finally, the MEM-prior is expected to yield more reliable 
results than the other type of entropy because our choice of the 
a priori image can closely describe a compact object. 

3. Description of the simulations 

We now describe all the various simulated data that we compute 
for different objects, (u, v) coverages, and signal-to-noise ratios 
(SNR), as well as the parameters used for the image reconstruc- 
tions. 

3.1. Simulated data 

Our simulated data sets are saved as OI-Fits files (Pauls et al. 
2005) and depend on several setting of the object type, the (u, v) 
coverage, and the SNR. The 90 data files are available on the 
JMMC website 1 . 

3.1 .1 . Astrophysical objects 

For our simulations, we consider ten astrophysical objects (see 
Fig. 1) that differ in term of their morphology and the typical 
length scales of their structures. 

1. LkHa: the model of LkHa describes a compact object with 
a peak of intensity and a smooth envelope. The model comes 
from the 2004 International Imaging Beauty Contest orga- 
nized by P. Lawson for the IAU (Lawson et al. 2004). 



1 http : //apps . jmmc . fr/oidata/shared/srenard/ 
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Fig. 1. Astrophysical objects used in the simulations. 
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2. A stellar surface: this is a model of the supergiant a Ori 
that presents some convective cells, producing small scales 
on a smooth background. This model was produced by A. 
Chiavassa for the 2009 Workshop on Interferometry Imaging 
(WII09) organized by J. -P. Berger and F. Malbet (Berger et 
al., in prep.). 

3. A stellar cluster: the model consists of a hundred point 
sources. The position and the brightness of the sources fol- 
low a normal law. 

4. Eta Carinae: the image of Eta Carinae presents many dif- 
ferent scales and structures, such as the extended gas and 
the stars. This image was retrieved from the Hubble Space 
Telescope's website 2 . Some treatments have been applied to 
the image, i.e. a mean of the three different color channels to 
produce a grayscale image and a cut of the low intensities to 
produce a null background. 

5. A protoplanetary disk: the model represents an 
HerbigAe/Be star with a point source (the star) and an 
extended structure (the disk). This model was computed by 
J. -P. Berger for the phase A science case of the VLTI-Spectro 
Imager instrument (Filho et al. 2008). 

6. A limb-darkened star: we used the power-law model of 
Hestroffer (1997) with an exponent a = 0.3. The image has 
a very smooth core with steep edges. 

7. The galaxy system M51: this image of M51 consists of 
as many different scales and structures as Eta Carinae (gas, 
stars, spiral arms). This image was retrieved from the Hubble 
Space Telescope's website and was processed in a way sim- 
ilar to the image of Eta Carinae. 

8. The AGN M87: this AGN has a jet that consists of a narrow 
structure surrounded by a smooth background due to the gas. 
The image was retrieved from the Hubble Space Telescope's 
website and the same treatments as the Eta Carinae 's image 
have been applied. 

9. A gravitational microlensing image: Gravitational mi- 
crolensing is an astronomical phenomenon due to the grav- 
itational lens effect. When a distant star or quasar becomes 
sufficiently aligned with a massive compact foreground ob- 
ject, the bending of light due to its gravitational field leads to 
two distorted unresolved images resulting in an observable 
magnification. The image shows four very compact struc- 
tures. This model was developed by J. Surdej for the phase A 
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Fig. 2. (u, v) coverage. From left to right: rich (245 sampled fre- 
quencies), medium (88 sampled frequencies), and poor (31 sam- 
pled frequencies). 



science case of the VLTI-Spectro Imager instrument (Filho 
et al. 2008). 

10. A rapid rotator: the rapid rotation of a star affects the stel- 
lar shape and the local emitted flux. We used the model of 
Domiciano de Souza et al. (2002), with parameters D = 0.78 
and T p = 35000 K. The resulting light distribution was pro- 
jected onto a plane with an inclination of 45°. 

Since the goal of our tests is to determine the influence of the 
object's type on the image reconstruction, all the considered ob- 
jects have a similar angular size of ~ 34 mas, which is consistent 
with the typical size of the field of view of an optical interferom- 
eter such as Amber/VLTI. To generalize our results, we expect 
the most important figure for a given object and instrument to be 
the number of resolved elements, which is equal to the ratio of 
the angular size of the object support to the effective resolution 
of the imaging system. Estimated by this ratio, the complexity of 
the objects that we have considered is in the range of 200 - 600 
resolved elements. 



3.1.2. (m,v) coverage 

To study the influence of the instrumental configuration and an- 
alyze the capability of the regularization to interpolate the avail- 
able data and thus fill the voids in the (w, v) plane, we consider 
several sparse (u, v) coverages (see Fig. 2). To remain general, 
our different (u, v) coverages were constructed to be uniformly 
spread. Each (u, v) coverage consists of concentric rings modu- 
lated by a sinusoid along the ring and with phases of the mod- 
ulations chosen to maximize the distance between the points of 
two adjacent rings. This also avoids symmetries in the (u, v) cov- 
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erage. The concentration of the rings is more important on small 
baselines than on the largest ones to insure a good sampling of 
low spatial frequencies. In this paper, we consider three (u, v) 
coverages which differ in terms of the number of samples, called 
hereafter rich (245 sampled frequencies), medium (88 sampled 
frequencies), and poor (31 sampled frequencies) coverages. The 
chosen (u, v) configurations could be considered as very good by 
the standards of existing data, but the goal of the paper is not to 
cover all possible cases but to show the general trend. 

3.1.3. Signal-to-noise ratio 

To investigate the effects of varying the SNR, we use a SNR 
factor y in the standard deviations given by 

o-g, m = j{Qm) and o~^ m = y . (14) 

Therefore with these settings, the error bars become 

0-//, m = 0~±,m =y(Qm). (15) 

To analyze the influence of the SNR on the reconstructed im- 
ages, three values of y are tested: high SNR (1%) , intermediate 
SNR (5%), and low SNR (10%). We limit our study to these SNR 
values to maintain a small amount of results. Therefore, there is 
no systematic attempt to search for the limit of SNR and we re- 
alize that our worst case (10%) might be considered moderate 
noise in real data. 

3.2. Parameters of the synthesized image 

The resolution and the field of view (hereafter FOV) of the re- 
constructed image have to be chosen with care. On the one 
hand, the pixel size A6 must be small enough to properly ac- 
count for the highest spatial frequencies available in the data. 
Using Shannon's rule, we find that A6 < A/(2B max ), where 
B max is the maximum length of the observed baselines and A 
the wavelength. In our reconstructions, the pixel size is one third 
of the limit set by Shannon's rule. This oversampling allows us 
to check whether image reconstruction can achieve some level 
of super-resolution. For instance, this yields AO = 0.4 mas at 
A = 2.2 \xm with B max - 190 m. On the other hand, the FOV has 
to be large enough to avoid field aliasing and therefore we take 
an image three times larger than the object itself. As all our ob- 
jects have a size of 34 mas and thus approximately fit into 85 x 85 
pixels, the reconstructed images have 256 x 256 pixels. 

3.3. Reconstruction strategy 

Given the data (determined by the object, the (u, v) coverage, 
and the noise realization) and the regularization, we conduct a 
sequence of image reconstructions for different values of the reg- 
ularization weight fi. Since the problems solved are convex (as 
explained in Sect. 2.3.1), their solutions do not depend on the 
starting point. To reduce the calculation time, we therefore try to 
use a starting point that is as close as possible to the final solu- 
tion. We begin the sequences of reconstructions with the highest 
value of fi and use the true image as the starting point for this first 
reconstruction. We then reduce \i and use the image previously 
obtained as the starting point. This latter step is repeated until 
we reach the lowest value of [i. Each reconstruction is an itera- 
tive process that is stopped when the norm of the gradient of the 
penalty function f(x) is below a preset threshold. This threshold 



is derived from the true image according to 
1/ (* ref )l 

|| W eC )|| < n = 10-5 \ U ( x ref)| (lfi) 

where n > is a small value and ^(jc re f) the penalty computed 
for the true image and for each value of [i. The £2 norm of the 
true image is normalized and we assume that n = 10" 5 . 

3.4. Image quality criterion: the mean-squared error 

To assess the quality of the reconstructed images, we consider 
the mean- squared error (hereafter MSE) of the reconstructed im- 
ages. In our simulations, we use normalized input images with 
the same pixel size as for the reconstructed image. Hence, to 
compare the reconstructed image x rec and the true image x ref , 
we can simply use the MSE defined as 

MSE = jf Z ( x ™ - = I n* rec - *ii 2 ' (17) 

n 

where TV is the total number of pixels. We note that, since we use 
complex visibilities and our priors, apart from the FOV one, are 
shift-invariant, the reconstructed image is correctly centered and 
there is no need to compensate for registration errors. 

4. Results and discussion 

We performed a total of -24000 simulations corresponding to 
the reconstruction of all cases described in Sect. 3 with different 
values of the hyperparameter fi. We present here the results of 
these simulations and analyze the consequences for the image 
reconstruction process in order to draw general conclusions. 

We begin by discussing the optimal value of the hyperpa- 
rameter and determining whether the MSE is a good quality cri- 
terion. We then discuss the effects of the following parameters: 

- the regularization: what are the good and bad regulariza- 
tions? 

- the limits: what are the minimal (w, v) coverage and SNR 
value? 

- the hyperparameter \i\ what is the optimal value? With which 
parameters does it vary? 

- the likelihood term: how does it vary? Can it be used to tune 
the regularization term instead of the hyperparameter fi! 

- the effective resolution: what degree of super-resolution can 
be achieved? How does it vary? 

4.1. Optimal regularization weight ji + 

We first investigate whether there is an optimal regularization 
weight \i for a given situation. Therefore, for each object, config- 
uration, SNR level, and regularization, we reconstruct an image 
for different values of the hyperparameter \i. 

The top row of the right panel of Fig. 3 shows the effects of 
different values of yu, whereas the left panel displays the obtained 
MSE (see Sect. 3.4) for each values of \i. For too small a value of 
yu, the under-regularized image (labeled with an A) has plenty of 
artifacts. In contrast, for too large a/i, the over-regularized image 
(labeled with a C) is blurred, and many fine features are lost. 
These two extreme situations correspond to high values of the 
MSE (A and C), but there is a situation where the MSE reaches 
a minimum and the image appears to have far fewer artifacts (B). 
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Fig. 3. Left panel shows a plot of the MSE as a function of the hyperparameter p. The different colors correspond to different 
levels of SNR (red high, blue intermediate, green poor). For each curve, the optimal value p + is labeled by a number (7, 2, and 3). 
The corresponding images are shown in the bottom row of the right panel. The top row of the right panel shows three reconstructed 
images with different values of p, labeled by a letter on the red curve of the left part (A an under-regularized image, B the best image, 
and C an over-regularized image). This example is made for the galaxy object and the medium (u, v) coverage. The regularization 
is the MEM-prior one. 



Visual comparison of the A and C images with the one obtained 
for B confirms that the MSE is a good criterion to correctly set 
the regularization weight. 

We conclude that there is indeed an optimal value of p, the 
one that gives the smallest MSE 



P 



argmin||x[f c -x ref || 2 



(18) 





where x|f c is the image reconstructed with a regularization 
weight set to p. As the minimum of the curve is quite flat, the 
optimal value of p is not precisely defined but may vary by a 
factor as large as either two or three with a negligible influence 
on the result. 

This procedure cannot be used in practice because the true 
image is obviously unknown. Nevertheless, this procedure al- 
lows us to define the most accurate image that can be recon- 
structed given the data and the type of regularization. In the 
following analysis, the reconstructed images are always ob- 
tained with p - p + . 

4.2. Dependence of p + on the MSE quality criterion 

We wish to determine the type of relationship that links the op- 
timized regularization factor p + to the MSE in order to detect 
any trends. The scatter plot in Fig. 4 reports the values of p + 
and MSE obtained for each simulation. In the left panel, the dif- 
ferent colors and symbols correspond to different objects. In the 
right panel, the different colors and symbols correspond to differ- 
ent regularizations. In the bottom row, we present our computed 
histograms of MSE (left) and p + (right). As in the rest of the pa- 
per, the plotted histograms are approximations of the probability 
density functions (PDF) of our results. These curves were com- 
puted from our samples following the optimal method described 
in Scott (1992). 

In the left part of Fig. 4, the different colors representing the 
objects seem to be aligned vertically. We therefore conclude that 
the MSE depends mostly on the structure of the observed object. 

More precisely, two classes of objects can be distinguished 
in the distribution of MSE in the left panel of Fig. 4. We there- 
fore grouped together the objects with similar behaviors in the 
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Fig. 5. Upper row: Distribution of the MSE (left) and the MSE + 
(right). The colors and letters represent the two classes of ob- 
jects: blue/B for the objects with very compact structures, red/C 
for the others. The total distribution is shown in the black/A 
curve. Bottom row: example of reconstructed images for the 
good (left) and bad (right) MSE + peak. 



top panels of Fig. 5: (i) the objects with very compact sources, 
i.e. the star cluster, the protoplanetary disk and the microlensing 
(curve in blue labeled B), and (ii) the other objects with extended 
structures (curve in red labeled C). The MSE is systematically 
higher for objects of the first class. 

Following these observations, we try to find a way of renor- 
malizing the MSE as a function of the object. To join the curves 
together, we define a new MSE criterion, called MSE + , by di- 
viding each MSE by the smallest MSE for each object sepa- 
rately. As expected, this normalization cancels the different ob- 
ject classes, as seen in the right part of Fig. 5. Two peaks clearly 
appear on the graph, distinguishing the good reconstructions 
(left peak) from the bad ones (right peak). An example of each 
case is shown on the bottom part of Fig. 5. The low quality re- 
constructions are caused by bad configurations (not enough (u, v) 
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Fig. 4. Top row: Scatter plots representing the optimal value of the hyperparameter [i + as a function of the MSE of the images. 
Bottom row: The corresponding histograms of MSE and Left column: The colors and symbols indicate the different classes of 
objects. Right column: The colors and symbols indicate the different regularization classes. 



points, low SNR) and/or bad regularizations as it can be seen in 
the right part of Fig. 6 and will be described in the next sections. 
To be sure that the peaks represent the good and bad reconstruc- 
tions and do not come from different object's types, we verify 
that each object is represented in each peak, as shown in the left 
part of Fig. 6. 

By visual inspection, we assess that the value of the MSE + 
leads to a correct ordering of the images reconstructed for a 
given object when the other settings change (data quality, type 
of regularization, etc.): the lower the MSE + , the higher the qual- 
ity of the image is. Other attempts at the renormalization of the 
MSE are explained in Appendix B. 

Now that a good quality criterion has been defined and that 
the optimal value of [i have been obtained, we can study the other 
parameters. 

4.3. Limits due to the (u, v) coverage and the SNR 

In this section, we classify the observational configurations 
((u, v) coverage and SNR) on the basis of the MSE + . For each 
data set (unique combination of object and regularization), we 
order the pair [(u, v) coverage, SNR] according to the value of 
MSE + they reach, giving them a rank from 1 for the best con- 
figuration (lowest MSE + ) to 9 for the worst (highest MSE + ). In 




MSE + MSE+ 

Fig. 6. Distribution of MSE + . Left: histograms of MSE + for dif- 
ferent objects in different colors; the gray zone corresponds to 
the total distribution, all objects confounded. Right: solid line, 
all the configurations and regularizations are kept; dashed line, 
with the sparsest (u, v) coverage removed; dot-dashed line, with 
the bad regularizations removed; in gray zone, with the sparsest 
(u, v) coverage and bad regularizations removed. 

the left panel of Fig. 7, we display the cumulative distributions 
of the ranks reached by every configuration, determining how 
many times a given configuration reaches at least the first rank, 
the second rank, etc. The highest quality configurations are the 
ones in the upper-left part of the plot. 
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Fig. 7. Left: cumulative distributions of the ranks reached by the different configurations of (w, v) coverage and SNR. Right: the 
histograms of the MSE + for different configurations of (w, v) coverage and SNR represented by different colors. 



The poor (u, v) coverage combined with any value of SNR is 
clearly too sparse to reconstruct good images. While acceptable 
reconstruction is possible for any considered SNR when there 
are enough samples in the (u, v) coverage. We deduce that there 
is a minimal (u, v) coverage needed to perform image recon- 
struction, whereas there is no such limit set by the tested SNR. 
However, when the (w, v) coverage is sufficient, the higher the 
SNR or the more filled the (u, v) coverage, the higher the quality 
of the reconstructed image is. The bottom row of the right panel 
of Fig. 3 shows how the visual quality of the optimal image de- 
pends on the SNR. As expected, the higher the SNR, the better 
the reconstructed image is. 

Figure 7 (right) shows the histogram of the MSE + for dif- 
ferent configurations of (u, v) coverage and SNR. It indicates 
that the success image reconstruction is more influenced by the 
amount of data than by the SNR: the MSE + is lower for a rich 
(u, v) coverage with a poor SNR than for a medium (u, v) cover- 
age with an intermediate SNR. We note that all the tested (u, v) 
coverage are uniform, but we expect that the amount of data has 
to be sufficient and also homogeneously spread in the (u, v) plan. 

After the removal of the sparsest coverage, the MSE + dis- 
tribution is shown in Fig. 6 in dashed line. There is still a little 
bump of bad MSE + caused by bad regularizations, as discussed 
in the next section. 

4.4. Quality of the regularizations 

Using the same principles as in the previous section, we clas- 
sify the regularizations as a function of the MSE + for different 
realizations. Figure 8 shows the corresponding cumulative dis- 
tributions. Isotropic TV appears to be the most successful regu- 
larization for the two classes of objects. The compactness prior 
with w[f or = \6 n \ 2 is the second highest quality choice for the 
very compact objects. The worst regularizations are the ^,-norm, 
the MEM-sqrt, and the MEM-log, as expected and explained in 
Sect. 2.3.2. Reconstructions for good and bad regularizations are 
illustrated in Fig. 9. 

In the MSE + distribution after the elimination of the bad reg- 
ularizations (see Fig. 6 in dot-dashed line), there is still much ev- 
idence of lower quality MSE + . We conclude that the bad MSE + 
are mainly caused by the sparsest (u, v) coverage. However, both 
have to be eliminated to obtain a cleaned sample of reconstructed 
images. 



In what follows, we made a selection and only retained 
the regularizations ( quadratic and l\ smoothness, compactness, 
isotropic TV and MEM with a Gaussian prior) and the (u, v) 
coverages (rich and medium) that lead to correct reconstructed 
images. We kept data for all values of SNR. 

4.5. Predetermined value of the hyperparameter ji + ? 

In a fully Bayesian framework, the data noise and the variabil- 
ity of the object are independent. We therefore expect that ///prior 
does not depend on the data (Tarantola 2005). As a result, for a 
given type of object, the optimal value for // should be the same 
regardless of the SNR or the (u, v) coverage. However, we are not 
in a truly Bayesian framework because the regularizations are 
derived from general considerations that must help us to solve 
for the degeneracies of the problem and are not really derived 
from the statistics of the brightness distributions of the observed 
objects. Since the degeneracies of the problem are mostly due 
to the sparsity of the (u, v) coverage, this observational parame- 
ter may have some influence on the tuning of the regularization 
weight. Our expectation is thus that, for a given type of object, 
a quasi optimal value of ji + can be derived from simulations and 
from simple considerations to scale this parameter if different 
image reconstruction settings are used. 

The hyperparameter ji + depends mostly on the regulariza- 
tions, as seen in the right part of Fig. 4, where the colors, repre- 
senting the regularizations, appear aligned horizontally and the 
pre-eminent peaks of the histograms are separated. Moreover, 
the hyperparameter // appears to be quasi independent of the 
SNR and the (u, v) coverage, as expected. As seen in Fig. 10 
for the TV regularization, if there is still a variation in the value 
of yu, it depends on the object morphology but not on the (w, v) 
plane and the SNR value. 

This finding allows us to link each regularization to a mean 
value of the hyperparameter //. This mean value gives a useful 
start point for the user (see Table 1). The equations to rescale 
this value in the case of different image settings are given in 
Appendix C. 
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Fig. 8. Cumulative distributions of the ranks reached by the regularizations. Left: all objects; Middle: objects with very small 
structures; Right: other objects. 
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regularizations. From left to right, TV, compactness 2 , £ p norm with p = 2, MEM-log. 
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Table 1. Table of the mean values of ji for each regularization. 



Regularizations 




Quadratic smoothness 


10* 


Compactness r 2 


10 7 


Compactness r 3 


10 7 


TV isotropic 


10 3 


i\ smoothness 


10 13 


MEM with prior 


10 2 



4.6. How different noise realizations affect the MSE and the 
optimal jd? 

In all simulations, in order not to influence the classification of 
the results, the same random seed is used to compute the noisy 



data. Therefore, we test 100 noise realizations for each regular- 
ization in the case of the galaxy with a medium (u, v) coverage 
and an intermediate SNR to study its impact on the curves com- 
puted from a single realization. From Fig. 11, we conclude that 
the MSE is not very different, thus the image reconstruction does 
not depend on the particular noise realization (as shown for ex- 
ample in Fig. 12). Moreover, the optimal value of \i varies by 
less than an order of magnitude. 

4. 7. The effective spatial resolution 

To investigate whether super-resolution is achievable and to 
quantify the amount of trustable super-resolution, we estimate 
the effective resolution of the reconstructed images. We define 
the effective resolution as the full width at half maximum (here- 
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Fig. 13. Left: FWHM of the Gaussian computed for the effective resolution in units of the interferometric resolution of the data. 
Up-right: comparison between the FWHM with (green) and without (magenta) the constraint of positivity. Bottom-right: variation 
in the effective resolution as a function of the hyperparameter ji for three different SNR (red high, blue intermediate, green poor). 
The regularization used is the TV one. 



A 
LU 
CO 



LU 2 
CO 



RGL quad, smooth 



RGL TV iso. 



RGL compact. 9 2 

i 


RGL compact. 3 


■ i _ 

RGL L1 smooth. 


RGL MEM prior 


I * 1 



10001 1.0 10001 

Hyperparameter li / < ji > 



Fig. 11. Variation in MSE and \i with different noise realizations. 
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50% in solid line, and 75% in dot line). Red, the mean optimal [i 
(cross) and its variance. 



after FWHM) of the Gaussian that yields the smallest value 
of MSE between the reconstructed image and the true image 
smoothed by that Gaussian 



FWHM = argmin ||x rec - G(w) • x 



refll 



(19) 



where G(w) is a linear smoothing filter that convolves its argu- 
ment by a Gaussian of FWHM equal to w. 

As shown in the left part of Fig. 13, the effective resolution 
varies with both the amount of data and the SNR: the more ex- 
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Fig. 12. Reconstructed images for the best (left) and the worst 
(right) noise realization in the TV case. 



tensively the (u, v) parameter space is filled and the higher the 
SNR, the higher the effective resolution is. As for the MSE (cf 
Sect. 4.3), the effective resolution is more influenced by the (u, v) 
coverage than by the SNR. Super-resolution can be achieved in 
only the best cases. 

The bottom-right part of Fig. 13 shows the expected behavior 
of the effective resolution: the higher the hyperparameter [i, the 
lower the effective resolution is. In order words, the more the 
image is regularized, the fewer the details that are visible. 

An interesting question is why this super-resolution exists. 
Although we do not have a definitive answer, we can speculate 
that it is due to the role of the positivity in the image reconstruc- 
tion (Biraud 1969). This is confirmed by our simulations (cf the 
upper right part of Fig. 13): without the positivity constraints, the 
distribution of the FWHMs has a peak higher than when using 
the positivity constraints. 
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Fig. 14. Histograms of //A/data an d /data/Afdata- •S'oZ/J /me: all configurations and regularizations are kept. Dashed line: with the 
sparsest (w, v) coverage removed. Dot-dashed line: with the bad regularizations removed. In gray zone: with the sparsest (w, v) 
coverage and bad regularizations removed. 



4.8. Other methods for tuning the regularization? 

For fully Gaussian statistics (i.e. all the likelihood and prior 
penalties are quadratic and there are no constraints such as the 
normalization and the positivity), the expected values of the to- 
tal penalty function f(x + ) and the likelihood term /data(* + ) at 
the optimal solution x + is given by (Tarantola 2005) 



E{f(x + )} = N 6a{a , 
E {/data(* + )} = Ndata " N e6i , 



(20) 
(21) 



where Afdata is the number of data points (for both visibility am- 
plitudes and phases), and Af e df is the number of equivalent de- 
grees of freedom. Despite our being unable to use fully Gaussian 
statistics, the left panel of Fig. 14 shows that the distribution of 
f(x + ) peaks approximatively at the value of A^ata- The spread 
of this distribution however prevents us to be able to tune the 
regularization level according to the criterion that f(x + ) « A^ata- 

The value of /data(* + )/Afdata is around 0.1, which means that 
A^edf /A^data ~ 90% of the data information is resolved. The im- 
age reconstruction is able to estimate almost as many parame- 
ters as data points. Since the /data(* + )/Ndata ratio has a smaller 
range (from -0.3 to -0.003) than the possible value of the hy- 
perparameter \i (from 100 to 10 13 ), it may be easier to tune the 
balance between the terms in the penalty function thanks to the 
/data(-^ + )/A^data value instead of [i. However, this criterion may 
be really variable with the noise statistics and a study of the 
variation in /data(* + )/Afdata with different noise statistics should 
be done before using it as a tune factor (Gull 1988; Pichon & 
Thiebaut 1998). 

Another way of determining the value of [i is the so-called 
L-curve. The L-curve is a log-log plot of the regularization term 
/prior versus the likelihood term /data for a range of the hyperpa- 
rameter ji. In the L-curve criterion, the regularization parameter 
fi is such that the corresponding point on the L-curve lies in the 
corner. This choice is motivated by the corner being the sepa- 
ration between the flat part where the solution is dominated by 



regularization errors and the vertical part where it is dominated 
by the perturbation errors (Hansen 2000). The correct behavior 
of the L-curve is confirmed by the simulations, as seen in Fig. 15 
(right): since the curves are plotted as a function of the highest 
quality image, they should cross in their corner, which is glob- 
ally the case. We note that the shape of the L-curve seems more 
complicated as there are at least two corners and not only one 
(see Fig. 15 left). The L-curve appears to be an appropriate tool 
for finding the optimal value of \i but a more general study has 
to be done to confirm its suitability (see comparison with GCV) 
and its practical implementation in an image reconstruction al- 
gorithm. 

5. Conclusions 

Thanks to the use of a flexible algorithm devoted to image re- 
construction in optical interferometry, we have performed a de- 
tailed study of the regularization terms. This study is the first one 
to compare such a number of regularizations on an equal foot- 
ing, i.e. with the same algorithm and using the same data type. 
Performing these systematic tests has allowed us to discuss the 
different parameters and terms in the image reconstruction and 
to extract some practical rules, which are summarized in the fol- 
lowing: 

1 . A minimal (u, v) coverage is necessary to reconstruct an im- 
age. Even if the image quality improves with the SNR, such 
a limit does not exist for the SNR. In other words, in the im- 
age reconstruction technique, increasing the number of tele- 
scopes is more interesting than constructing larger ones. The 
homogeneity of the (w, v) coverage is probably also critical 
but has not been tested. 

2. Some regularizations are suitable for the optical image re- 
construction and others not, regardless of the object being 
targeted. The holes in the (u, v) plane are the major issue in 
optical interferometry and the main role of the regularization 
is the correct interpolation of the missing data, regardless 
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of the object. The highest quality regularization among the 
tested ones is the isotropic total variation. 

3. The hyperparameter \i does not depend on the (u, v) cov- 
erage and the SNR as theoretically expected and depends 
mostly on the type of regularization. An optimal value for 
each regularization tested in this paper is given in Table 1. 
This optimal value may vary by a factor of 2-3 without there 
being any major changes in the images. A slight dependence 
on both the object structures and pixel size is also discernible 
and the equation to rescale the optimal values are computed. 
It should be interesting to implement regularizations inde- 
pendent of the pixel size. 

4. Super-resolution can be achieved in the image reconstruc- 
tion process and its level rises with the (u, v) coverage filling. 

5. There are various possible ways of tuning the regularization 
level: 

- The visual tuning is enough as the \i value can slightly 
vary without causing any large changes in the recon- 
structed image. 

- Setting the likelihood term /data seems to be a more effec- 
tive way of fixing the balance between the regularization 
and the likelihood terms. However, the variation in the 
likelihood term with noise statistics needs to be investi- 
gated. 

- The L-curve criterion could give correct results. 
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Appendix A: The regularization expressions 

In our simulations, we test 1 1 different regularization terms com- 
monly used in image reconstruction methods and implemented 
inMiRA: 



1. Quadratic smoothness: 

/prior(-*0 — II-* — S • 



(A.1) 
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where S is a smoothing operator implemented via finite dif- 
ferences. 



2-3. Compactness (Le Besnerais et al. 2008): 

/prior W = J] W n m " X n ' 



(A.2) 



which is a separable quadratic regularization. To enforce 
compactness, the weights w{f' or > have to increase with 
the distance from the center of the image. Two different 
cases were studied in the simulations: w{f' or = \\0 n /A6\\ 2 and 
= \\0 n /A6\\ 3 , where ||0 n /A0|| is the distance in pixels of 



w 



prior _ 



the n-th pixel from the center of the field of view. 
4. Total variation (Strong & Chan 2003): 



/prior (X) = J] ^WVXn^lf + e 2 , (A3) 



where 



is the squared magnitude of the spatial gradient in the image, 
6 > is a small number inserted to avoid the discontinuity in 
zero, and (n\,ni) ~ n are the two dimensional indices of the 
n-th pixel. In our reconstructions, e has always been chosen 
so as to be negligible compared to the gradient of significant 
structures. 

5. — smoothness (Charbonnier et al. 1997): 

/priori) = T 2 J] <A(II( D • *)J/T) , (A.4) 

n 

where if/(z) = z - log(l + z) is a norm, r > is a thresh- 
old level, and D is a finite difference operator approximating 
the g-th spatial derivatives of its argument. When ||(D • x) n || 
is much smaller than the threshold, t 2 if/(\\(D • x) n \\/r) « 
1 /2 ||(D • x) n \\ 2 , while, for ||(D • it is much larger than the 
threshold, r 2 ^(||(D • x) n \\/r) « t||(D • x) n \\. This regulariza- 
tion attempts to strongly smooth the weak spatial gradients 
and slightly smooth the strong gradients. In our simulations, 
we take D to approximate the spatial Laplacian and choose a 
threshold small enough such that the regularization behaved 
mostly like a linear smoothness. 



6-8. ^-norm: 



/prlorW = 2 n (^ + e2 r /2 *Zn 



(A.5) 



where e > is a small value introduced to avoid the singular- 
ity in zero when p < 1. For p > 1, the ^,-norm regularization 
tends to produce a smooth image as it reduces the variance 
of the pixels. For p < 1, the ^-norm regularization tends 
to promote sparsity in the image. This is interesting mostly 
for objects consisting of point sources. True sparsity con- 
straints would be obtained for p = 0, although when p < 1 
the regularization is no longer convex and the optimization 
problem becomes extremely difficult to solve as p gets closer 
to 0. Results in compress sensing (Candes et al. 2006) have 
proven that choosing p - 1 yields the most sparse solution, 
like p - for a large class of problems. However, taking 
p - 1 yields a non-smooth but convex problem that is much 
easier to solve than the combinatorial problem resulting from 





1st normalized MSE 



2nd normalized MSE 



Fig. B.l. Attempt to renormalize the MSE. Distribution of the 
first normalized MSE (left) and of the second normalized MSE 
(right). The colors and letters represent the two classes of ob- 
jects: blue/B for the objects with very compact structures, red/C 
for the others. The total distribution is shown in black/A curve. 

the choice p = 0. With positivity and normalization con- 
straints, the ^i-norm of x is constant. Hence, taking p = 1 is 
meaningless in our framework and we consider only p - 1.5, 
p - 2 and p - 3. 

9-11. Maximum entropy methods (Narayan & Nityananda 
1986): 

/prior (X) = ~ h ( X n '> X n) • ( A -6) 

Here the prior is to assume that the image is drawn towards a 
prior model x according to a non-quadratic potential h, called 
the entropy. We try three entropies: 

MEM-sqrt: h(x\ x) = V* ; (A.7) 
MEM-log: h(x; x) = log(x) ; (A.8) 
MEM-prior: h(x; x) = x - x - x log (x/x) . (A.9) 

MEM was first introduced in radioastronomy and is useful 
for images made of bright point-like sources on a smooth 
background. In our simulations, we took the prior image x 
to be the isotropic 2-D Gaussian that most accurately fits the 
amplitude visibility data. 

Appendix B: Renormalization of MSE 

We attempted to normalize the MSE with two additional meth- 
ods: 

1. Since the squared difference between the real image and a 
smoothed version of the real image is higher for images with 
sharp or point-like structures, we computed a first normal- 
ized MSE as 



MSE 



norm.,1 



2j« (x n x n ) 



(B.l) 



where S is a smoothing operator. The distribution of this nor- 
malized MSE is shown in the right panel of Fig. B. 1 : the dis- 
tribution is narrower but the two classes remain, despite the 
normalization. 

In the second normalization, the MSE is compared to the 
norm of the reference image 



MSE 



norm., 2 



2j« (*n X n ) 



(B.2) 
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This normalization is more effective than the previous one, 
as it joins the curves together. However, the distribution is 
unimodal and does not enable us to distinguish the good and 
the bad reconstructions. It is thus not a useful normalization. 



Appendix C: Scaling the hyperparameter \i 

In a Bayesian framework, the prior penalty [i / pr j 0r (x) should only 
depend on both the sought after brightness distribution 1(0) and 
the image parametrization. In this appendix and from these sim- 
ple principles, we derive a method to adapt the value of the hy- 
perparameter ji when the image parameters such as the pixel size 
are modified. 

Using the sampled image model in Eq. (4), the normalization 
of x implies that 

where we have used the Riemann approximation of integrals, AO 
is the pixel size, and L = 2 is the number of dimensions in 0. 
Without any loss of generality, we can assume that the sought 
after brightness distribution is normalized such that J 1(0) d L = 
1 ; hence 

a*(A0) L , (C.l) 



and Eq. (4) becomes 



C.1. Separable £ p norm 



(A6) L I(0 n ). 



(C.2) 



Using Eq. (4) and then, the Riemann approximation, the prior 
penalty for a separable £ p norm regularization, after substituting 
in Eq. (A. 5), becomes 

V /prior (X) =^1^ ^ *^ P Yj |/( ^ )|/? 



in 



d0, 



(C.3) 



which shows that for a regularization by a separable £ p norm 

(A0 L /a p in general; ^ 
li OC { ' _ . b (C.4) 
[A0 {p) with the normalization constraint. 

Hence, with the normalization constraint, the optimal value of \i 
should be the same regardless of the pixel size for a regulariza- 
tion given by a separable £\ norm. 



C.2. l p norm on the gradient 

Using 1-D notation to simplify the equations, the prior penalty 
for a regularization by the £ p norm on the gradient is given by 

A*/prior(*) = V \ Xn+l ~ Xn ^ 

n 

xfia>> YjW n + Ae)-I(6 n )\P 

n 

^ fia p A0P YjW(O n )\ p , 



where del (6) denotes the partial derivative of the brightness dis- 
tribution along the angular direction. In L dimensions and using 
the Riemann approximation, this gives 

n /prior (x)^ l ia p A6 p - L f \d G I(0)\ p dO, 



which shows that 



A6 L p /a p in general; 
^ y/\ i L -p( L+l ) with the normalization constraint. 

Applying this result for a regularization by quadratic smooth- 
ness in 2D, e.g. Eq. (A.l), we found that, with a normalization 
constraint, \i oc A0~ 4 . 



C.3. Total variation 

The preceding result, with p = 1, readily applies to regulariza- 
tion by the total variation, that is 



A6 L 1 la in general; 



AO 



i-i 



with the normalization constraint. 



(C.6) 



We can also deduce that, if a relaxed version of TV is used, as 
in Eq. (A.3), the relaxation parameter e must scale as the pixel 
size AO to have the prior penalty approximately insensitive to the 
pixel size. 

We note that, with our particular choice of the threshold r 
for the £2 - £\ smoothness regularization defined in Eq. (A.4), 
we expect this regularization to behave mostly like TV. 

C.4. Quadratic compactness 

The quadratic compactness we used in MiRA is given by 
Eq. (A.2) 



li /prior (X) = n Yj X ' * * Z ^ ^n) 



\A6 



2 

AW 



AQq+L 



l(0f dO, 



with q = 2 or 3. From this last approximation, we derive the 
scaling of [i with the pixel size 



jd oc 



A6 q /a in general; 



A6 q ~ 



with the normalization constraint. 



(C.7) 



Hence, in 2D (L = 2) and for a normalized image, [i does not 
depend on the pixel size for q = 2, while it scales as A# for 
q = 3. 

C.5. Maximum entropy 
For maximum entropy methods, we have 
_ 0,1/2 r 

^^^n^ K0) ll2 d0, 

n 

and 

!i Yj h ^n\ x n ) * fi J h(l(0); 1(0)) d0 , 
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with h(x; x) = x-x-x log(x/x) as in Eq. (A. 9). from which we 
deduce that 

( A6 L /a l/2 for MEM-sqrt; 
^ 00 \A6 L /a for MEM-prior. ( ' } 

Hence, for a normalized image, [i does not depend on the pixel 
size in a MEM-prior regularization. 



